GLMM per art
Torgeir
18 februar, 2021
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.3 v purrr 0.3.4
## v tibble 3.0.6 v dplyr 1.0.4
## v tidyr 1.1.2 v stringr 1.4.0
## v readr 1.4.0 v forcats 0.5.0## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(ggeffects) # Estimated Marginal Means and Marginal Effects from Regression Models
# more at: https://strengejacke.github.io/ggeffects/
library(performance) # diagnostic-plots to check assumptions
library(report) # Result-summaries in text-format## report is in alpha - help us improve by reporting bugs on github.com/easystats/report/issues# Data drom Data_exploration2_nesting.R
time.dep <- readRDS("time.dep") %>%
mutate(time.deploy = time.deploy/10) # shortening / zooming out on deploy
# to avoid warnings during model fittingPurpose
This notebook is meant to set up the final GLMModels going into my thesis. Upset is identical to the “glmm_in_process”-file, but structure is different as I’m going to make one model for each species.
Formula
The model formula I will use is $n time.deploy * flash $ for each species, and my \(\alpha = 0.05\).
sp <- c("raadyr", "rev", "hjort", "grevling", "elg", "gaupe")
time.dep2c <- time.dep %>%
rename(species = validated_species) %>% #shortening name
filter(species %in% sp) %>% #filtering out species
# including Control as part of the flash-column, since it differs from flash=0
mutate(flash = factor(
ifelse(period == "Control", "Control", flash)),
week = lubridate::isoweek(date))
time.dep2 <- time.dep2c %>%
filter(!period == "Control") %>% # but removing it for now, because it is set up in a longer timeframe
mutate(flash = factor(flash, levels = c(0,1)) )
class(time.dep2$flash)## [1] "factor"Not all periods have identical length. Hence, I need to set a maximum length for my period durations. As proposed by Neri, I will calculate the median for white LED-periods and IR-periods, and use the smallest median to shorten all periods overextending that value.
First I’ll filter out any periods shorter than 4 days ( as of 18.02.2021, only 1 period ). Then I’ll cut the duration of all periods overextending the smallest median.
# filter out shortest periods, and find median period length
cut <- 0.4 # setting 4 days as the minimum length of a period
# find lengths
time.period <- time.dep2 %>% group_by(loc, period, flash) %>%
summarise(period_length = max(time.deploy))
# checking which periods will be removed
time.period %>% filter(period_length < cut) %>%
arrange(period_length) #%>% kableExtra::kable("html")## # A tibble: 1 x 4
## # Groups: loc, period [1]
## loc period flash period_length
## <fct> <chr> <fct> <dbl>
## 1 829 1_1 1 0# then merge lengths and filter based on that
time.dep3 <- time.dep2 %>% left_join(time.period) %>%
filter(period_length > cut)
# find median length after filtering short lengths out
time.period %>% filter(period_length > cut) %>% filter(flash == 1) %>%
summary() # median period length 85 days, mean: 84## loc period flash period_length
## 15 : 2 Length:62 0: 0 Min. : 0.800
## 127 : 2 Class :character 1:62 1st Qu.: 6.600
## 193 : 2 Mode :character Median : 8.500
## 231 : 2 Mean : 8.421
## 257 : 2 3rd Qu.:10.850
## 456 : 2 Max. :13.200
## (Other):50time.period %>% filter(period_length > cut) %>% filter(flash == 0) %>%
summary() # median period length 79 days, mean: 89## loc period flash period_length
## 15 : 2 Length:63 0:63 Min. : 1.200
## 127 : 2 Class :character 1: 0 1st Qu.: 6.550
## 193 : 2 Mode :character Median : 7.900
## 231 : 2 Mean : 8.946
## 257 : 2 3rd Qu.:11.850
## 456 : 2 Max. :18.300
## (Other):51# extract lengths of each unique period
h <- time.dep3 %>% group_by(loc, period, period_length, flash)%>% nest() %>%
select(!data)
#extracting median and multiplying by 10, to use in the correctly scaled plot
hh <- h$period_length[h$flash == 1] %>% median() # median white LED
hh <- c(hh, h$period_length[h$flash == 0] %>% median()) * 10 # + median IR
# smallest median
h <- min(hh)79 days is the median period length of IR-periods, and it is shorter than the median of white LED flash. The summary also tells us that
Periods below 4 days has now been removed, which as of 18.02.2021 (before updated data from Neri) were 1 period, which was a flash period of 0 days: |loc| period|flash |period_length| |—|——-|——|:———–:| |829| 1_1 | 1 | 0.0 |
# plot periods with median as intercept
p_td <- time.dep3 %>%
ggplot(aes(loc, 10*time.deploy, colour = period, )) +
geom_line(aes(linetype = flash),position = position_dodge(width = 1), lineend = "square") +
coord_flip() +
labs(title = "Period lengths per camera",
x = "Location", y = "Time since deployment",
caption = "Dotted lines reprecent median period length for IR and white LED.\n Data superceding that is trimmed away for the GLMM-modelling.") +
ggpubr::theme_classic2() #+ theme(legend.position = "right") find way to set legend inside
p_td + geom_hline(aes(yintercept = h), linetype = "dashed", alpha =.5) +
geom_hline(aes(yintercept = max(hh)), linetype = "dashed", alpha =.5) +
#annotate(geom = "text",x=4, y=h+8.6, label = "- median", size = 3, alpha =.7) +
scale_y_continuous(breaks = sort(c(0, 50, hh, 100, 150)))
# failed attempts that could inspire a better plot later
# p_td + geom_hline(aes(yintercept = h))+ # using median days as intercept
# scale_y_continuous(breaks = sort(c(ggplot_build(p_td)$layout$panel_ranges[[1]]$y.major_source, h)))
# geom_text(aes(25, h, label = "median", vjust = -1), nudge_y = 10, show.legend = F)There was an overweight of IR-periods extending past the median line.
Lastly, performing the filter:
# filtering out periods longer than (shortest) median length.
time.dep4 <- time.dep3 %>% filter(time.deploy < h/10)Modelling
Roe deer
Diagnostics
# filter species
sp = "raadyr"
time_sp <- filter(time.dep4, species %in% sp) #.dep4 = trimmed data
# Model
m_sp <- glmer(n.obs ~ time.deploy * flash + # fixed effects
(1 | loc) + (1 | week), # random effects
data = time_sp,
family = poisson) # poisson family of distributions
# ggpredict is similar to expand.grid
p_sp <- ggeffects::ggpredict(m_sp, terms = c("time.deploy", "flash"))
# Diagnostics
plot(p_sp, add.data = TRUE) + labs(caption = "add.data = TRUE")
plot(p_sp, residuals = TRUE) + labs(caption = "residuals")
performance::check_model(m_sp) # check assumptions## Loading required namespace: qqplotr## `geom_smooth()` using formula 'y ~ x'## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 8696 rows containing missing values (geom_text_repel).## `geom_smooth()` using formula 'y ~ x'## `geom_smooth()` using formula 'y ~ x'## Warning: ggrepel: 128 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
The response variable is a summary of number of events per day. Most days had no roe deer. In the performance-test for model assumptions it is clear that some assumptions aren’t met.
Negative
In the non-normality of resiudals-plot, the residuals skew off from the line when moving towards positive quantiles. There is not a homogeneity of variance. Maybe this could be fixed by centering the n.obs-column?
There appeares to be five influential observations in the Cook’s distance-plot, maybe more, as the warning from ggrepel refers to 93 unlabeled data points due to overlaps.
Positive
Although the model has an interaction term between flash and time since deployment, the multicollinearity between them is low! My random effects follows a normal distribution.
Concluding
I am not sure about how to make up for breaking these assumptions. For now, I will go on completing models for the rest of the species.
# Summary, report, model
summary(m_sp)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: n.obs ~ time.deploy * flash + (1 | loc) + (1 | week)
## Data: time_sp
##
## AIC BIC logLik deviance df.resid
## 5273.9 5316.5 -2631.0 5261.9 8823
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3051 -0.3109 -0.1816 -0.0961 13.4925
##
## Random effects:
## Groups Name Variance Std.Dev.
## week (Intercept) 0.3147 0.561
## loc (Intercept) 2.2252 1.492
## Number of obs: 8829, groups: week, 52; loc, 33
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.349052 0.299975 -11.164 <2e-16 ***
## time.deploy 0.008207 0.025277 0.325 0.745
## flash1 0.169515 0.126476 1.340 0.180
## time.deploy:flash1 -0.035958 0.029867 -1.204 0.229
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tm.dpl flash1
## time.deploy -0.304
## flash1 -0.228 0.565
## tm.dply:fl1 0.199 -0.663 -0.842report::report(m_sp) # text-summary of my model, to include in a report## We fitted a poisson mixed model (estimated using ML and Nelder-Mead optimizer) to predict n.obs with time.deploy and flash (formula: n.obs ~ time.deploy * flash). The model included loc and week as random effects (formula: list(~1 | loc, ~1 | week)). The model's total explanatory power is substantial (conditional R2 = 0.43) and the part related to the fixed effects alone (marginal R2) is of 4.11e-04. The model's intercept, corresponding to time.deploy = 0 and flash = 0, is at -3.35 (95% CI [-3.94, -2.76], p < .001). Within this model:
##
## - The effect of time.deploy is non-significantly positive (beta = 8.21e-03, 95% CI [-0.04, 0.06], p = 0.745; Std. beta = 0.02, 95% CI [-0.09, 0.13])
## - The effect of flash [1] is non-significantly positive (beta = 0.17, 95% CI [-0.08, 0.42], p = 0.180; Std. beta = 0.04, 95% CI [-0.10, 0.17])
## - The interaction effect of flash [1] on time.deploy is non-significantly negative (beta = -0.04, 95% CI [-0.09, 0.02], p = 0.229; Std. beta = -0.08, 95% CI [-0.21, 0.05])
##
## Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using the Wald approximation.plot(p_sp)
Model interpretation
The intercept-value is considered significantly negative, which is to say that there were a low chance of detecting any roe deer at an IR-camera the same day I visited the camera.
I saw a roe deer about to walk by a CT when I came to inspect it. The roe deer saw me and fled, right before it was detected by the camera. Chances are I’ve scared animals other times as well, but haven’t noticed it.
The effect of Time since deployment is non-significant, and \(\beta = 0.008\). That means there is no difference on the baseline detection rate for an IR camera over time (after controlling for seasonal changes).
For white LED flash \(\beta = 0.170\), meaning that the intercept is slightly higher than for IR, but the difference is non-significant (\(p = 0.18\)). However, the detection rate is slightly decreasing the longer the white LED stays, which differs from the IR, but the effect is non-significant (\(p = 0.23\)).
Hypothesising
If there truly is an effect of the white LED for long periods on the detection rate of roe deer, this effect could in turn account for the different intercept values of IR and flash, as the IR periods often start after a flash period.
Remembering my study design, 20 cameras start with white LED, 20 with IR. Intercept should be equal (1st period). 2nd period; white LED moved, new LED CTs (same intercept), new IR CTs (hypothetical lower intercept due to flash effect). 3rd period; white LED moved, new LED CTs (IR intercept), new IR CTs (hypothetical lower intercept). 4th period; - - || - - , new LED CTs (- - | | - - ), new IR CTs ( - - | | - - ).
Which sums up to 3 IR periods where detection rates could start lower than that of white LED. And over time the lack of white LED flash in the new IR sites would account for an increase in detection rates.
The effect would of course vary because some locations experienced gaps due to full SD cards or empty batteries.
f
Red Fox
# filter species
sp = "rev"
time_sp <- filter(time.dep4, species %in% sp) #.dep4 = trimmed data
# Model
m_sp <- glmer(n.obs ~ time.deploy * flash + # fixed effects
(1 | loc) + (1 | week), # random effects
data = time_sp,
family = poisson) # poisson family of distributions
# ggpredict is similar to expand.grid
p_sp <- ggeffects::ggpredict(m_sp, terms = c("time.deploy", "flash"))
# Diagnostics
plot(p_sp, add.data = TRUE) + labs(caption = "add.data = TRUE")
plot(p_sp, residuals = TRUE) + labs(caption = "residuals")
performance::check_model(m_sp) # check assumptions## `geom_smooth()` using formula 'y ~ x'## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 8627 rows containing missing values (geom_text_repel).## `geom_smooth()` using formula 'y ~ x'## `geom_smooth()` using formula 'y ~ x'## Warning: ggrepel: 198 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Interpret
# Summary, report, model
summary(m_sp)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: n.obs ~ time.deploy * flash + (1 | loc) + (1 | week)
## Data: time_sp
##
## AIC BIC logLik deviance df.resid
## 3264.1 3306.7 -1626.1 3252.1 8823
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.7743 -0.2355 -0.1734 -0.1301 12.2970
##
## Random effects:
## Groups Name Variance Std.Dev.
## week (Intercept) 0.1284 0.3583
## loc (Intercept) 0.6869 0.8288
## Number of obs: 8829, groups: week, 52; loc, 33
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.53418 0.21574 -16.381 <2e-16 ***
## time.deploy 0.01304 0.03447 0.378 0.705
## flash1 0.13452 0.18606 0.723 0.470
## time.deploy:flash1 -0.01509 0.04274 -0.353 0.724
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tm.dpl flash1
## time.deploy -0.585
## flash1 -0.460 0.567
## tm.dply:fl1 0.400 -0.663 -0.854report::report(m_sp) # text-summary of my model, to include in a report## We fitted a poisson mixed model (estimated using ML and Nelder-Mead optimizer) to predict n.obs with time.deploy and flash (formula: n.obs ~ time.deploy * flash). The model included loc and week as random effects (formula: list(~1 | loc, ~1 | week)). The model's total explanatory power is moderate (conditional R2 = 0.19) and the part related to the fixed effects alone (marginal R2) is of 4.68e-04. The model's intercept, corresponding to time.deploy = 0 and flash = 0, is at -3.53 (95% CI [-3.96, -3.11], p < .001). Within this model:
##
## - The effect of time.deploy is non-significantly positive (beta = 0.01, 95% CI [-0.05, 0.08], p = 0.705; Std. beta = 0.03, 95% CI [-0.12, 0.18])
## - The effect of flash [1] is non-significantly positive (beta = 0.13, 95% CI [-0.23, 0.50], p = 0.470; Std. beta = 0.08, 95% CI [-0.11, 0.27])
## - The interaction effect of flash [1] on time.deploy is non-significantly negative (beta = -0.02, 95% CI [-0.10, 0.07], p = 0.724; Std. beta = -0.03, 95% CI [-0.22, 0.15])
##
## Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using the Wald approximation.plot(p_sp)
Badger
# filter species
sp = "grevling"
time_sp <- filter(time.dep4, species %in% sp) #.dep4 = trimmed data
# Model
m_sp <- glmer(n.obs ~ time.deploy * flash + # fixed effects
(1 | loc) + (1 | week), # random effects
data = time_sp,
family = poisson) # poisson family of distributions
# ggpredict is similar to expand.grid
p_sp <- ggeffects::ggpredict(m_sp, terms = c("time.deploy", "flash"))
# Diagnostics
plot(p_sp, add.data = TRUE) + labs(caption = "add.data = TRUE")
plot(p_sp, residuals = TRUE) + labs(caption = "residuals")
performance::check_model(m_sp) # check assumptions## `geom_smooth()` using formula 'y ~ x'## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 8561 rows containing missing values (geom_text_repel).## `geom_smooth()` using formula 'y ~ x'## `geom_smooth()` using formula 'y ~ x'## Warning: ggrepel: 267 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Interpret
# Summary, report, model
summary(m_sp)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: n.obs ~ time.deploy * flash + (1 | loc) + (1 | week)
## Data: time_sp
##
## AIC BIC logLik deviance df.resid
## 3064.1 3106.7 -1526.1 3052.1 8823
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.9115 -0.2138 -0.1318 -0.0726 14.5242
##
## Random effects:
## Groups Name Variance Std.Dev.
## week (Intercept) 1.601 1.265
## loc (Intercept) 1.180 1.086
## Number of obs: 8829, groups: week, 52; loc, 33
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.391474 0.320432 -13.705 <2e-16 ***
## time.deploy 0.069170 0.037446 1.847 0.0647 .
## flash1 -0.061728 0.197962 -0.312 0.7552
## time.deploy:flash1 0.007123 0.043589 0.163 0.8702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tm.dpl flash1
## time.deploy -0.453
## flash1 -0.336 0.551
## tm.dply:fl1 0.316 -0.651 -0.875report::report(m_sp) # text-summary of my model, to include in a report## We fitted a poisson mixed model (estimated using ML and Nelder-Mead optimizer) to predict n.obs with time.deploy and flash (formula: n.obs ~ time.deploy * flash). The model included loc and week as random effects (formula: list(~1 | loc, ~1 | week)). The model's total explanatory power is substantial (conditional R2 = 0.40) and the part related to the fixed effects alone (marginal R2) is of 3.77e-03. The model's intercept, corresponding to time.deploy = 0 and flash = 0, is at -4.39 (95% CI [-5.02, -3.76], p < .001). Within this model:
##
## - The effect of time.deploy is non-significantly positive (beta = 0.07, 95% CI [-4.22e-03, 0.14], p = 0.065; Std. beta = 0.15, 95% CI [-9.36e-03, 0.32])
## - The effect of flash [1] is non-significantly negative (beta = -0.06, 95% CI [-0.45, 0.33], p = 0.755; Std. beta = -0.04, 95% CI [-0.23, 0.15])
## - The interaction effect of flash [1] on time.deploy is non-significantly positive (beta = 7.12e-03, 95% CI [-0.08, 0.09], p = 0.870; Std. beta = 0.02, 95% CI [-0.17, 0.21])
##
## Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using the Wald approximation.plot(p_sp)
Moose
# filter species
sp = "elg"
time_sp <- filter(time.dep4, species %in% sp) #.dep4 = trimmed data
# Model
m_sp <- glmer(n.obs ~ time.deploy * flash + # fixed effects
(1 | loc) + (1 | week), # random effects
data = time_sp,
family = poisson) # poisson family of distributions## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0842096 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?# ggpredict is similar to expand.grid
p_sp <- ggeffects::ggpredict(m_sp, terms = c("time.deploy", "flash"))
# Diagnostics
plot(p_sp, add.data = TRUE) + labs(caption = "add.data = TRUE")
plot(p_sp, residuals = TRUE) + labs(caption = "residuals")
performance::check_model(m_sp) # check assumptions## `geom_smooth()` using formula 'y ~ x'## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 8505 rows containing missing values (geom_text_repel).## `geom_smooth()` using formula 'y ~ x'## `geom_smooth()` using formula 'y ~ x'## Warning: ggrepel: 323 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Interpret
# Summary, report, model
summary(m_sp)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: n.obs ~ time.deploy * flash + (1 | loc) + (1 | week)
## Data: time_sp
##
## AIC BIC logLik deviance df.resid
## 2049.6 2092.1 -1018.8 2037.6 8823
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.5446 -0.1784 -0.1215 -0.0680 18.3565
##
## Random effects:
## Groups Name Variance Std.Dev.
## week (Intercept) 0.4087 0.6393
## loc (Intercept) 1.6108 1.2692
## Number of obs: 8829, groups: week, 52; loc, 33
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.5119563 0.0008386 -5380.386 < 2e-16 ***
## time.deploy -0.0066476 0.0008394 -7.919 2.39e-15 ***
## flash1 0.3354212 0.0008385 400.045 < 2e-16 ***
## time.deploy:flash1 -0.0449704 0.0008388 -53.610 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tm.dpl flash1
## time.deploy 0.000
## flash1 0.000 0.000
## tm.dply:fl1 0.000 0.001 0.000
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0842096 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?report::report(m_sp) # text-summary of my model, to include in a report## We fitted a poisson mixed model (estimated using ML and Nelder-Mead optimizer) to predict n.obs with time.deploy and flash (formula: n.obs ~ time.deploy * flash). The model included loc and week as random effects (formula: list(~1 | loc, ~1 | week)). The model's total explanatory power is substantial (conditional R2 = 0.31) and the part related to the fixed effects alone (marginal R2) is of 2.15e-03. The model's intercept, corresponding to time.deploy = 0 and flash = 0, is at -4.51 (95% CI [-4.51, -4.51], p < .001). Within this model:
##
## - The effect of time.deploy is significantly negative (beta = -6.65e-03, 95% CI [-8.29e-03, -5.00e-03], p < .001; Std. beta = -0.01, 95% CI [-0.23, 0.20])
## - The effect of flash [1] is significantly positive (beta = 0.34, 95% CI [0.33, 0.34], p < .001; Std. beta = 0.17, 95% CI [-0.08, 0.42])
## - The interaction effect of flash [1] on time.deploy is significantly negative (beta = -0.04, 95% CI [-0.05, -0.04], p < .001; Std. beta = -0.10, 95% CI [-0.35, 0.16])
##
## Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using the Wald approximation.plot(p_sp)
Red deer
# filter species
sp = "hjort"
time_sp <- filter(time.dep4, species %in% sp) #.dep4 = trimmed data
# Model
m_sp <- glmer(n.obs ~ time.deploy * flash + # fixed effects
(1 | loc) + (1 | week), # random effects
data = time_sp,
family = poisson) # poisson family of distributions
# ggpredict is similar to expand.grid
p_sp <- ggeffects::ggpredict(m_sp, terms = c("time.deploy", "flash"))
# Diagnostics
plot(p_sp, add.data = TRUE) + labs(caption = "add.data = TRUE")
plot(p_sp, residuals = TRUE) + labs(caption = "residuals")
performance::check_model(m_sp) # check assumptions## `geom_smooth()` using formula 'y ~ x'## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 8701 rows containing missing values (geom_text_repel).## `geom_smooth()` using formula 'y ~ x'## `geom_smooth()` using formula 'y ~ x'## Warning: ggrepel: 127 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Interpret
# Summary, report, model
summary(m_sp)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: n.obs ~ time.deploy * flash + (1 | loc) + (1 | week)
## Data: time_sp
##
## AIC BIC logLik deviance df.resid
## 1153.7 1196.2 -570.8 1141.7 8823
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.6201 -0.1121 -0.0399 -0.0252 24.4856
##
## Random effects:
## Groups Name Variance Std.Dev.
## week (Intercept) 0.4017 0.6338
## loc (Intercept) 4.7039 2.1688
## Number of obs: 8829, groups: week, 52; loc, 33
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.87201 0.63691 -9.220 < 2e-16 ***
## time.deploy -0.14656 0.08114 -1.806 0.07088 .
## flash1 -0.78251 0.43274 -1.808 0.07057 .
## time.deploy:flash1 0.31983 0.09856 3.245 0.00117 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tm.dpl flash1
## time.deploy -0.381
## flash1 -0.318 0.496
## tm.dply:fl1 0.298 -0.715 -0.858report::report(m_sp) # text-summary of my model, to include in a report## We fitted a poisson mixed model (estimated using ML and Nelder-Mead optimizer) to predict n.obs with time.deploy and flash (formula: n.obs ~ time.deploy * flash). The model included loc and week as random effects (formula: list(~1 | loc, ~1 | week)). The model's total explanatory power is substantial (conditional R2 = 0.46) and the part related to the fixed effects alone (marginal R2) is of 0.01. The model's intercept, corresponding to time.deploy = 0 and flash = 0, is at -5.87 (95% CI [-7.12, -4.62], p < .001). Within this model:
##
## - The effect of time.deploy is non-significantly negative (beta = -0.15, 95% CI [-0.31, 0.01], p = 0.071; Std. beta = -0.32, 95% CI [-0.68, 0.03])
## - The effect of flash [1] is non-significantly negative (beta = -0.78, 95% CI [-1.63, 0.07], p = 0.071; Std. beta = 0.38, 95% CI [-0.05, 0.82])
## - The interaction effect of flash [1] on time.deploy is significantly positive (beta = 0.32, 95% CI [0.13, 0.51], p < .01; Std. beta = 0.71, 95% CI [0.28, 1.14])
##
## Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using the Wald approximation.plot(p_sp)
Lynx
# filter species
sp = "gaupe"
time_sp <- filter(time.dep4, species %in% sp) #.dep4 = trimmed data
# Model
m_sp <- glmer(n.obs ~ time.deploy * flash + # fixed effects
(1 | loc) + (1 | week), # random effects
data = time_sp,
family = poisson) # poisson family of distributions
# ggpredict is similar to expand.grid
p_sp <- ggeffects::ggpredict(m_sp, terms = c("time.deploy", "flash"))
# Diagnostics
plot(p_sp, add.data = TRUE) + labs(caption = "add.data = TRUE")
plot(p_sp, residuals = TRUE) + labs(caption = "residuals")
performance::check_model(m_sp) # check assumptions## `geom_smooth()` using formula 'y ~ x'## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 8306 rows containing missing values (geom_text_repel).## `geom_smooth()` using formula 'y ~ x'## `geom_smooth()` using formula 'y ~ x'## Warning: ggrepel: 521 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Interpret
# Summary, report, model
summary(m_sp)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: n.obs ~ time.deploy * flash + (1 | loc) + (1 | week)
## Data: time_sp
##
## AIC BIC logLik deviance df.resid
## 483.8 526.3 -235.9 471.8 8823
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.3592 -0.0579 -0.0355 -0.0246 25.3606
##
## Random effects:
## Groups Name Variance Std.Dev.
## week (Intercept) 0.9442 0.9717
## loc (Intercept) 1.9729 1.4046
## Number of obs: 8829, groups: week, 52; loc, 33
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.7403 0.7999 -9.676 <2e-16 ***
## time.deploy 0.1663 0.1286 1.294 0.1958
## flash1 1.4764 0.7035 2.099 0.0358 *
## time.deploy:flash1 -0.2160 0.1535 -1.407 0.1593
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tm.dpl flash1
## time.deploy -0.729
## flash1 -0.646 0.719
## tm.dply:fl1 0.528 -0.750 -0.862report::report(m_sp) # text-summary of my model, to include in a report## We fitted a poisson mixed model (estimated using ML and Nelder-Mead optimizer) to predict n.obs with time.deploy and flash (formula: n.obs ~ time.deploy * flash). The model included loc and week as random effects (formula: list(~1 | loc, ~1 | week)). The model's total explanatory power is substantial (conditional R2 = 0.32) and the part related to the fixed effects alone (marginal R2) is of 0.02. The model's intercept, corresponding to time.deploy = 0 and flash = 0, is at -7.74 (95% CI [-9.31, -6.17], p < .001). Within this model:
##
## - The effect of time.deploy is non-significantly positive (beta = 0.17, 95% CI [-0.09, 0.42], p = 0.196; Std. beta = 0.37, 95% CI [-0.19, 0.93])
## - The effect of flash [1] is significantly positive (beta = 1.48, 95% CI [0.10, 2.86], p < .05; Std. beta = 0.69, 95% CI [-0.02, 1.39])
## - The interaction effect of flash [1] on time.deploy is non-significantly negative (beta = -0.22, 95% CI [-0.52, 0.08], p = 0.159; Std. beta = -0.48, 95% CI [-1.15, 0.19])
##
## Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using the Wald approximation.plot(p_sp)
Hare, deer and squirrelywhere (and pine marten)
Blueprint for other species in chunk below:
Diagnostics
# filter species
sp = "xx"
time_sp <- filter(time.dep4, species %in% sp) #.dep4 = trimmed data
# Model
m_sp <- glmer(n.obs ~ time.deploy * flash + # fixed effects
(1 | loc) + (1 | week), # random effects
data = time_sp,
family = poisson) # poisson family of distributions
# ggpredict is similar to expand.grid
p_sp <- ggeffects::ggpredict(m_sp, terms = c("time.deploy", "flash"))
# Diagnostics
plot(p_sp, add.data = TRUE) + labs(caption = "add.data = TRUE")
plot(p_sp, residuals = TRUE) + labs(caption = "residuals")
performance::check_model(m_sp) # check assumptionsInterpret
# Summary, report, model
summary(m_sp)
report::report(m_sp) # text-summary of my model, to include in a report
plot(p_sp)